# SETUP ----

#setting working directory- CHANGE TO YOUR WORKING DIRECTORY
setwd()

#loading necessary packages
library(conflicted)
library(dplyr)
conflicts_prefer(dplyr::summarize)
conflicts_prefer(dplyr::filter)
library(ggplot2)
library(ggpubr)
library(graphics)
library(haven)
library(lmtest)
library(Matrix)
library(PanelMatch)
library(panelView)
library(plm)
library(readr)
library(readxl)
library(sandwich)
library(stargazer)
library(tidyverse)
library(texreg)

# IMPORTING AND CLEANING DATA ----

#Loading data for insurgent conscription
dat <- read_csv('ICdata_jcr.csv')

#removing X1 column
IC = subset(dat, select = c(2:34))

names(IC)

lapply(IC, class)

#making sure variables are of the correct class
IC$dyadid <- as.numeric(IC$dyadid)
IC$dyadid <- as.integer(IC$dyadid)
IC$year <-as.integer(IC$year)
IC$trt <- as.numeric(as.character(IC$trt))
IC$alttrt <- as.numeric(as.character(IC$alttrt))
IC$icabandon <- as.numeric(IC$icabandon)
IC$icadopt<- as.numeric(IC$icadopt)
IC<- as.data.frame(IC)

#Loading data for analysis with RHRV
rhrv_cont <- read_csv('rhrv_cont.csv')
rhrv_cont <- rhrv_cont[,c(2:71)]

rhrv_cont$dyadid <- as.integer(rhrv_cont$dyadid)
rhrv_cont$year <-as.integer(rhrv_cont$year)
rhrv_cont$trt <- as.numeric(as.character(rhrv_cont$trt))
rhrv_cont$excluded <- as.numeric(as.character(rhrv_cont$excluded))
rhrv_cont$exclpopn <- as.numeric(as.character(rhrv_cont$exclpopn))
rhrv_cont$dyadsfirstyearconflict <- as.numeric(as.character(rhrv_cont$dyadsfirstyearconflict))
rhrv_cont$fcap <- as.numeric(as.character(rhrv_cont$fcap))
rhrv_cont$rsup <- as.numeric(as.character(rhrv_cont$rsup))
rhrv_cont$irregular2 <- as.numeric(as.character(rhrv_cont$irregular2))
rhrv_cont$intensity <- as.numeric(as.character(rhrv_cont$intensity))
rhrv_cont$ideolcomm <- as.numeric(as.character(rhrv_cont$ideolcomm))

rhrv_cont <- as.data.frame(rhrv_cont)

# REPLICATION DATA FOR MAIN RESULTS----

## REPLICATION FOR TABLE 2: Frequency of Insurgent Conscription  ----

# Calculating descriptives 

# 1) Total number of unique armed groups
total_unique_armed_groups <- length(unique(IC$dyadid))

# 2) Count of unique armed groups where ic is ever and never equal to 1
armed_groups_ever_ic_1 <- IC %>%
  group_by(dyadid) %>%
  summarize(ic_ever = any(ic == 1), ic_never = all(ic != 1)) %>%
  summarize(
    total_unique_dyads = n (),
    ever_equal_1 = sum(ic_ever),
    never_equal_1 = sum(ic_never)
  )


# 3) Percentage of unique armed groups where ic is ever and never equal to 1
percent_armed_groups <- armed_groups_ever_ic_1 %>%
  mutate(
    percent_ever = ever_equal_1 / total_unique_armed_groups * 100,
    percent_never = never_equal_1 / total_unique_armed_groups * 100
  )

# 4) Total number of unique years
total_unique_years <- length(unique(IC$year))

# 5) Count of unique years where ic is ever and never equal to 1
years_ever_ic_1 <- IC %>%
  group_by(year) %>%
  summarize(ic_ever = max(ic == 1), ic_never =  all(ic != 1)) %>%
  summarize(
    total_unique_years = n (),
    ever_equal_1 = sum(ic_ever),
    never_equal_1 = sum(ic_never)
  )

# 6) Percentage of unique years where ic is ever and never equal to 1
percent_years <- years_ever_ic_1 %>%
  mutate(
    percent_ever = ever_equal_1 / total_unique_years * 100,
    percent_never = never_equal_1 / total_unique_years * 100
  )

# 7) Total number of unique countries
total_unique_countries <- length(unique(IC$ccode))

# 8) Count of unique countries where ic is ever and never equal to 1
countries_ever_ic_1 <- IC %>%
  group_by(ccode) %>%
  summarize(ic_ever = max(ic == 1), ic_never =  all(ic != 1)) %>%
  summarize(
    total_unique_countries = n (),
    ever_equal_1 = sum(ic_ever),
    never_equal_1 = sum(ic_never)
  )

# 9) Percentage of unique countries where ic is ever and never equal to 1
percent_countries <- countries_ever_ic_1 %>%
  mutate(
    percent_ever = ever_equal_1 / total_unique_countries * 100,
    percent_never = never_equal_1 / total_unique_countries * 100
  )

# 10) Total number of unique dyadid-year combinations
total_unique_dyadid_years <- nrow(unique(IC[, c("dyadid", "year")]))

# Count of unique dyadid-year combinations where ic is 1 and 0
dyadid_years_ic <- IC %>%
  group_by(dyadid, year) %>%
  summarize(ic_ever = any(ic == 1), ic_never = all(ic != 1)) %>%
  ungroup() %>%
  summarize(
    total_dyadid_years = n(),
    ever_equal_1 = sum(ic_ever),
    never_equal_1 = sum(ic_never)
  )

# Percentage of unique dyadid-year combinations where ic is 1 and 0
percent_dyadid_years <- dyadid_years_ic %>%
  mutate(
    percent_ever = ever_equal_1 / total_dyadid_years * 100,
    percent_never = never_equal_1 / total_dyadid_years * 100
  )

#Inputted this info into overleaf table

## REPLICATION FOR TABLE 3: ATT of Sustained State Violence on Insurgent Conscription by Period ----

##PANEL 1

# These models estimate the effect of sustained targeting on insurgent conscription at time t+1- t+3 with 1 year lag and matching on missing. 
set.seed(12345)
m1 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 1, forbid.treatment.reversal = FALSE)
mset.m1 <- m1$att
PE.m1  <- PanelEstimate(sets = m1, data = IC)


set.seed(12345)
m2 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 2, forbid.treatment.reversal = FALSE)
mset.m2 <- m2$att
PE.m2  <- PanelEstimate(sets = m2, data = IC)

set.seed(12345)
m3 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 3, forbid.treatment.reversal = FALSE)

mset.m3  <- m3$att
PE.m3   <- PanelEstimate(sets = m3, data = IC)

summary(PE.m1)
summary(PE.m2)
summary(PE.m3)

##PANEL 2

# These models estimate the effect of sustained targeting on insurgent conscription at time t+1- t+3 with 1 year lag and no matching on missing. 
set.seed(12345)
m4 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = FALSE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 1, forbid.treatment.reversal = FALSE)
mset.m4 <- m4$att
PE.m4  <- PanelEstimate(sets = m4, data = IC)

set.seed(12345)
m5 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = FALSE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +   I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 2, forbid.treatment.reversal = FALSE)
mset.m5 <- m5$att
PE.m5  <- PanelEstimate(sets = m5, data = IC)

set.seed(12345)
m6 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = FALSE,
                 covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                 +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                 size.match = 15, qoi = "att"
                 , outcome.var = "ic",
                 lead = 3, forbid.treatment.reversal = FALSE)

mset.m6  <- m6$att
PE.m6   <- PanelEstimate(sets = m6, data = IC)

summary(PE.m4)
summary(PE.m5)
summary(PE.m6)

##PANEL 3
# These models estimate the effect of sustained targeting on insurgent conscription at time t+1- t+3 with 2 year lag and matching on missing. 
set.seed(12345)
m7 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                 +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic, 1:2)),
                 size.match = 10, qoi = "att"
                 , outcome.var = "ic",
                 lead = 1, forbid.treatment.reversal = FALSE)

mset.m7 <- m7$att
PE.m7 <- PanelEstimate(sets = m7, data = IC)

set.seed(12345)
m8 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                 +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic, 1:2)),
                 size.match = 10, qoi = "att"
                 , outcome.var = "ic",
                 lead = 2, forbid.treatment.reversal = FALSE)

mset.m8 <- m8$att
PE.m8 <- PanelEstimate(sets = m8, data = IC)

set.seed(12345)
m9 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                 treatment = "trt", refinement.method = "CBPS.weight",
                 data = IC, match.missing = TRUE,
                 covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                 +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic, 1:2)),
                 qoi = "att"
                 , outcome.var = "ic",
                 lead = 3, forbid.treatment.reversal = FALSE)

mset.m9 <- m9$att
PE.m9 <- PanelEstimate(sets = m9, data = IC)

summary(PE.m7)
summary(PE.m8)
summary(PE.m9)


## REPLICATION FOR FIGURE 1 ----
#This plots the results of Panel 1 of Table 3

png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR/icplot.png")
par(mfrow=c(2,2))
plot(PE.m1, ylab = "Estimated Effect",
     xlab = "Time", main = NA, ylim=c(0.00,0.30))
plot(PE.m2, ylab = "Estimated Effect",
     xlab = "Time", main = NA, ylim=c(0.00,0.30))
plot(PE.m3, ylab = "Estimated Effect",
     xlab = "Time", main = NA, ylim=c(0.00,0.30))
dev.off()

###REPLICATION FOR TABLE 4: Linear Probability Models

### creating variable for pre 1980

IC <- IC %>% mutate(pre1980 = ifelse(year<=1979, 1, 0))

###regression models

summary(lm_ic <- lm(ic  ~ trt, data=IC))
summary(lm_ic_sc <- lm(ic  ~ trt + dyadsfirstyearconflict + intensitynewnum + pre1980 + irregular2 +ideolcomm, data=IC))
summary(lm_ic_ac <- lm(ic ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 + ideolcomm +intensitynewnum + pre1980 + nat_resource_funding, data=IC))


###regression clustered on dyad
ic_d <- coeftest(lm_ic, vcov = vcovCL, cluster = ~dyadid)
ic_d_sc <- coeftest(lm_ic_sc, vcov = vcovCL, cluster = ~dyadid)
ic_d_ac <- coeftest(lm_ic_ac, vcov = vcovCL, cluster = ~dyadid)



###regression clustered on country

#ic 
ic_c <- coeftest(lm_ic, vcov = vcovCL, cluster = ~ccode)
ic_c_sc <- coeftest(lm_ic_sc, vcov = vcovCL, cluster = ~ccode)
ic_c_ac <- coeftest(lm_ic_ac, vcov = vcovCL, cluster = ~ccode)


##Table 4
stargazer(ic_d, ic_d_sc, ic_d_ac, ic_c, ic_c_sc, ic_c_ac, keep.stat=c("n"), 
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Pre-1980", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))


###REPLICATION FOR TABLE 5: State violence across levels of intensity

intensitytable <- IC %>%
  group_by(intensitynewnum) %>%
  summarize(
    count = n(),
    trt1_count = sum(trt == 1),
    trt1_percent = (sum(trt == 1) / n()) * 100
  )


# REPLICATION FOR APPENDICES ----

## APPENDIX B1: Clearest cases of IC ----

#Create indicator for clear case
IC <- IC %>% mutate(ic_clear = ifelse(ic==1 & difficult_case==0, 1, 0))

### Table B1 ----

##PANEL 1
set.seed(12345)
m1_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 1, forbid.treatment.reversal = FALSE)
mset.m1_icp <- m1_icp$att
PE.m1_icp  <- PanelEstimate(sets = m1_icp, data = IC)

set.seed(12345)
m2_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 2, forbid.treatment.reversal = FALSE)
mset.m2_icp <- m2_icp$att
PE.m2_icp  <- PanelEstimate(sets = m2_icp, data = IC)

set.seed(12345)
m3_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 3, forbid.treatment.reversal = FALSE)

mset.m3_icp  <- m3_icp$att
PE.m3_icp   <- PanelEstimate(sets = m3_icp, data = IC)

summary(PE.m1_icp)
summary(PE.m2_icp)
summary(PE.m3_icp)

##PANEL 2

set.seed(12345)
m4_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = FALSE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 1, forbid.treatment.reversal = FALSE)
mset.m4_icp <- m4_icp$att
PE.m4_icp  <- PanelEstimate(sets = m4_icp, data = IC)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+2 with 1 year lag and no matching on missing. It generates the results for row 2 of panel 1 in Table 2.
set.seed(12345)
m5_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = FALSE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +   I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 2, forbid.treatment.reversal = FALSE)
mset.m5_icp <- m5_icp$att
PE.m5_icp  <- PanelEstimate(sets = m5_icp, data = IC)

set.seed(12345)
m6_icp <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = FALSE,
                     covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                     +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic_clear, 1:1)),
                     size.match = 15, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 3, forbid.treatment.reversal = FALSE)

mset.m6_icp  <- m6_icp$att
PE.m6_icp   <- PanelEstimate(sets = m6_icp, data = IC)

summary(PE.m4_icp)
summary(PE.m5_icp)
summary(PE.m6_icp)

##PANEL 3
set.seed(12345)
m7_icp <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                     +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic_clear, 1:2)),
                     size.match = 10, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 1, forbid.treatment.reversal = FALSE)

mset.m7_icp <- m7_icp$att
PE.m7_icp <- PanelEstimate(sets = m7_icp, data = IC)

set.seed(12345)
m8_icp <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                     +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic_clear, 1:2)),
                     size.match = 10, qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 2, forbid.treatment.reversal = FALSE)

mset.m8_icp <- m8_icp$att
PE.m8_icp <- PanelEstimate(sets = m8_icp, data = IC)

set.seed(12345)
m9_icp <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                     treatment = "trt", refinement.method = "CBPS.weight",
                     data = IC, match.missing = TRUE,
                     covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                     +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(ic_clear, 1:2)),
                     qoi = "att"
                     , outcome.var = "ic_clear",
                     lead = 3, forbid.treatment.reversal = FALSE)

mset.m9_icp <- m9_icp$att
PE.m9_icp <- PanelEstimate(sets = m9_icp, data = IC)

summary(PE.m7_icp)
summary(PE.m8_icp)
summary(PE.m9_icp)

### Table B2 ----
###regression models

summary(lm_icp <- lm(ic_clear ~ trt, data=IC))
summary(lm_icp_sc <- lm(ic_clear  ~ trt + dyadsfirstyearconflict + intensitynewnum + pre1980 + irregular2 +ideolcomm, data=IC))
summary(lm_icp_ac <- lm(ic_clear ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + pre1980 + nat_resource_funding, data=IC))


###regression clustered on dyad
icp_d <- coeftest(lm_icp, vcov = vcovCL, cluster = ~dyadid)
icp_d_sc <- coeftest(lm_icp_sc, vcov = vcovCL, cluster = ~dyadid)
icp_d_ac <- coeftest(lm_icp_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
icp_c <- coeftest(lm_icp, vcov = vcovCL, cluster = ~ccode)
icp_c_sc <- coeftest(lm_icp_sc, vcov = vcovCL, cluster = ~ccode)
icp_c_ac <- coeftest(lm_icp_ac, vcov = vcovCL, cluster = ~ccode)


## Table B2


stargazer(icp_d, icp_d_sc, icp_d_ac, icp_c, icp_c_sc, icp_c_ac, keep.stat=c("n"), 
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Pre-1980", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))

## APPENDIX B2: including rebel goal in controls ----


### Table B3: ATT of Sustained State Violence on Insurgent Conscription by Period, with rebel goals ----

##PANEL 1
set.seed(12345)
PM.sec1 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~ I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) +  I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                      
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 1, forbid.treatment.reversal = FALSE)

PE.sec1  <- PanelEstimate(sets = PM.sec1, data = IC)

set.seed(12345)
PM.sec2 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(sec, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                      
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 2, forbid.treatment.reversal = FALSE)

PE.sec2  <- PanelEstimate(sets = PM.sec2, data = IC)

set.seed(12345)
PM.sec3 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(sec, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                      
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 3, forbid.treatment.reversal = FALSE)


PE.sec3  <- PanelEstimate(sets = PM.sec3, data = IC)

summary(PE.sec1)
summary(PE.sec2)
summary(PE.sec3)

##PANEL 2

set.seed(12345)
PM.sec1_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                         treatment = "trt", refinement.method = "CBPS.weight",
                         data = IC, match.missing = FALSE,
                         covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) +  I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(sec, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                         
                         size.match = 15, qoi = "att"
                         , outcome.var = "ic",
                         lead = 1, forbid.treatment.reversal = FALSE)


PE.sec1_nm  <- PanelEstimate(sets = PM.sec1_nm, data = IC)


set.seed(12345)
PM.sec2_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                         treatment = "trt", refinement.method = "CBPS.weight",
                         data = IC, match.missing = FALSE,
                         covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) +  I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(sec, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                         
                         size.match = 15, qoi = "att"
                         , outcome.var = "ic",
                         lead = 2, forbid.treatment.reversal = FALSE)


PE.sec2_nm  <- PanelEstimate(sets = PM.sec2_nm, data = IC)

set.seed(12345)
PM.sec3_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                         treatment = "trt", refinement.method = "CBPS.weight",
                         data = IC, match.missing = FALSE,
                         covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(sec, 1:1)) + I(lag(intensitynewnum, 1:1))+ I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                         
                         size.match = 15, qoi = "att"
                         , outcome.var = "ic",
                         lead = 3, forbid.treatment.reversal = FALSE)


PE.sec3_nm  <- PanelEstimate(sets = PM.sec3_nm, data = IC)

summary(PE.sec1_nm)
summary(PE.sec2_nm)
summary(PE.sec3_nm)

##PANEL 3

set.seed(12345)
PM.sec2_1 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                        treatment = "trt", refinement.method = "CBPS.weight",
                        data = IC, match.missing = TRUE,
                        covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) +  I(lag(fcap, 1:2)) +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(sec, 1:2)) + I(lag(intensitynewnum, 1:2))+ I(lag(nat_resource_funding, 1:2)) +  I(lag(ic, 1:2)),
                        
                        size.match = 10, qoi = "att"
                        , outcome.var = "ic",
                        lead = 1, forbid.treatment.reversal = FALSE)


PE.sec2_1  <- PanelEstimate(sets = PM.sec2_1 , data = IC)

set.seed(12345)
PM.sec2_2 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                        treatment = "trt", refinement.method = "CBPS.weight",
                        data = IC, match.missing = TRUE,
                        covs.formula = ~ I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) +  I(lag(fcap, 1:2)) +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(sec, 1:2)) + I(lag(intensitynewnum, 1:2))+  I(lag(nat_resource_funding, 1:2)) +  I(lag(ic, 1:2)),
                        size.match = 10, qoi = "att"
                        , outcome.var = "ic",
                        lead = 2, forbid.treatment.reversal = FALSE)

PE.sec2_2  <- PanelEstimate(sets = PM.sec2_2 , data = IC)

set.seed(12345)
PM.sec2_3 <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                        treatment = "trt", refinement.method = "CBPS.weight",
                        data = IC, match.missing = TRUE,
                        covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) +  I(lag(fcap, 1:2)) +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(sec, 1:2)) + I(lag(intensitynewnum, 1:2))+  I(lag(nat_resource_funding, 1:2)) +  I(lag(ic, 1:2)),
                        size.match = 10, qoi = "att"
                        , outcome.var = "ic",
                        lead = 3, forbid.treatment.reversal = FALSE)


PE.sec2_3  <- PanelEstimate(sets = PM.sec2_3 , data = IC)

summary(PE.sec2_1)
summary(PE.sec2_2)
summary(PE.sec2_3)

### Table B4 ----

###regression clustered on dyad

#ic
summary(lm_sec_sc <- lm(ic ~ trt + dyadsfirstyearconflict + intensitynewnum + pre1980 + irregular2 +ideolcomm + sec, data=IC))
summary(lm_sec_ac <- lm(ic ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + pre1980 + nat_resource_funding + sec, data=IC))


###regression clustered on dyad
sec_sc_d <- coeftest(lm_sec_sc, vcov = vcovCL, cluster = ~dyadid)
sec_ac_d <- coeftest(lm_sec_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
sec_sc_c<- coeftest(lm_sec_sc, vcov = vcovCL, cluster = ~ccode)
sec_ac_c <- coeftest(lm_sec_ac, vcov = vcovCL, cluster = ~ccode)

# table
stargazer(sec_sc_d, sec_ac_d, sec_sc_c, sec_ac_c, keep.stat=c("n"),
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Pre-1980", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Secessionist", "Constant"))
## APPENDIX B3: Replication for alternative treatment threshold ----

### Table B5 ----
##PANEL 1

set.seed(12345)
am1 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = TRUE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 1, forbid.treatment.reversal = FALSE)
mset.am1 <- am1$att
PE.am1  <- PanelEstimate(sets = am1, data = IC)

set.seed(12345)
am2 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = TRUE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 2, forbid.treatment.reversal = FALSE)
mset.am2 <- am2$att
PE.am2  <- PanelEstimate(sets = am2, data = IC)


set.seed(12345)
am3 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = TRUE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 3, forbid.treatment.reversal = FALSE)

mset.am3  <- am3$att
PE.am3   <- PanelEstimate(sets = am3, data = IC)

summary(PE.am1)
summary(PE.am2)
summary(PE.am3)

##PANEL 2


set.seed(12345)
am4 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = FALSE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 1, forbid.treatment.reversal = FALSE)
mset.am4 <- am4$att
PE.am4  <- PanelEstimate(sets = am4, data = IC)


set.seed(12345)
am5 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = FALSE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 2, forbid.treatment.reversal = FALSE)
mset.am5 <- am5$att
PE.am5  <- PanelEstimate(sets = am5, data = IC)

set.seed(12345)
am6 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "alttrt", refinement.method = "CBPS.weight",
                  data = IC, match.missing = FALSE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(ic, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "ic",
                  lead = 3, forbid.treatment.reversal = FALSE)

mset.am6  <- am6$att
PE.am6   <- PanelEstimate(sets = am6, data = IC)

summary(PE.am4)
summary(PE.am5)
summary(PE.am6)

### Table B6 ----

summary(lm_alt <- lm(ic ~ alttrt, data=IC))
summary(lm_alt_sc <- lm(ic  ~ alttrt + dyadsfirstyearconflict + intensitynewnum + pre1980 + irregular2 +ideolcomm, data=IC))
summary(lm_alt_ac <- lm(ic ~ alttrt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + pre1980 + nat_resource_funding, data=IC))


###regression clustered on dyad
alt_d <- coeftest(lm_alt, vcov = vcovCL, cluster = ~dyadid)
alt_d_sc <- coeftest(lm_alt_sc, vcov = vcovCL, cluster = ~dyadid)
alt_d_ac <- coeftest(lm_alt_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
alt_c <- coeftest(lm_alt, vcov = vcovCL, cluster = ~ccode)
alt_c_sc <- coeftest(lm_alt_sc, vcov = vcovCL, cluster = ~ccode)
alt_c_ac <- coeftest(lm_alt_ac, vcov = vcovCL, cluster = ~ccode)


##replication for Table 4


stargazer(alt_d, alt_d_sc, alt_d_ac, alt_c, alt_c_sc, alt_c_ac, keep.stat=c("n"), 
          covariate.labels=c("1 year state targeting threshold", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Pre-1980", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))

## APPENDIX B4: Replication with insurgent conscription adoption ----

### Table B7 ----

#PANEL 1
set.seed(12345)
PM.t1 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                    treatment = "trt", refinement.method = "CBPS.weight",
                    data = IC, match.missing = TRUE,
                    covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                    size.match = 10, qoi = "att"
                    , outcome.var = "icadopt",
                    lead = 1, forbid.treatment.reversal = FALSE)
mset.t1 <- PM.t1$att
PE.t1  <- PanelEstimate(sets = PM.t1, data = IC)

set.seed(12345)
PM.t1_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                       treatment = "trt", refinement.method = "CBPS.weight",
                       data = IC, match.missing = FALSE,
                       covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                       size.match = 10, qoi = "att"
                       , outcome.var = "icadopt",
                       lead = 1, forbid.treatment.reversal = FALSE)


mset.t1_nm <- PM.t1_nm$att
PE.t1_nm  <- PanelEstimate(sets = PM.t1_nm, data = IC)

summary(PE.t1)
summary(PE.t1_nm)

# not in analysis- just curious if rebels adopt at different durations following state violence

set.seed(12345)
PM.t2 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                    treatment = "trt", refinement.method = "CBPS.weight",
                    data = IC, match.missing = TRUE,
                    covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                    size.match = 10, qoi = "att"
                    , outcome.var = "icadopt",
                    lead = 2, forbid.treatment.reversal = FALSE)
mset.t2 <- PM.t2$att
PE.t2  <- PanelEstimate(sets = PM.t2, data = IC)

set.seed(12345)
PM.t2_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                       treatment = "trt", refinement.method = "CBPS.weight",
                       data = IC, match.missing = FALSE,
                       covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                       size.match = 10, qoi = "att"
                       , outcome.var = "icadopt",
                       lead = 2, forbid.treatment.reversal = FALSE)


mset.t2_nm <- PM.t2_nm$att
PE.t2_nm  <- PanelEstimate(sets = PM.t2_nm, data = IC)

summary(PE.t2)
summary(PE.t2_nm)

set.seed(12345)
PM.t3 <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                    treatment = "trt", refinement.method = "CBPS.weight",
                    data = IC, match.missing = TRUE,
                    covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                    size.match = 10, qoi = "att"
                    , outcome.var = "icadopt",
                    lead = 3, forbid.treatment.reversal = FALSE)
mset.t3 <- PM.t3$att
PE.t3  <- PanelEstimate(sets = PM.t3, data = IC)

set.seed(12345)
PM.t3_nm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                       treatment = "trt", refinement.method = "CBPS.weight",
                       data = IC, match.missing = FALSE,
                       covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1)),
                       size.match = 10, qoi = "att"
                       , outcome.var = "icadopt",
                       lead = 3, forbid.treatment.reversal = FALSE)


mset.t3_nm <- PM.t3_nm$att
PE.t3_nm  <- PanelEstimate(sets = PM.t3_nm, data = IC)

summary(PE.t3)
summary(PE.t3_nm)
### Table B8 ----

summary(lm_adopt <- lm(icadopt ~ trt, data=IC))
summary(lm_adopt_sc <- lm(icadopt  ~trt + dyadsfirstyearconflict + intensitynewnum + pre1980 + irregular2 +ideolcomm, data=IC))
summary(lm_adopt_ac <- lm(icadopt ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + pre1980 + nat_resource_funding, data=IC))


###regression clustered on dyad
adopt_d <- coeftest(lm_adopt, vcov = vcovCL, cluster = ~dyadid)
adopt_d_sc <- coeftest(lm_adopt_sc, vcov = vcovCL, cluster = ~dyadid)
adopt_d_ac <- coeftest(lm_adopt_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
adopt_c <- coeftest(lm_adopt, vcov = vcovCL, cluster = ~ccode)
adopt_c_sc <- coeftest(lm_adopt_sc, vcov = vcovCL, cluster = ~ccode)
adopt_c_ac <- coeftest(lm_adopt_ac, vcov = vcovCL, cluster = ~ccode)

stargazer(adopt_d, adopt_d_sc, adopt_d_ac, adopt_c, adopt_c_sc, adopt_c_ac, keep.stat=c("n"), 
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Pre-1980", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))

## APPENDIX C: Replication for RHRV Models ----

### Table C2 ----

set.seed(12345)
m1_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 1, forbid.treatment.reversal = FALSE)
mset.m1_rhrv <- m1_rhrv$att
PE.m1_rhrv  <- PanelEstimate(sets = m1_rhrv, data = rhrv_cont)



#This model estimates the effect of sustained targeting on insurgent conscription at time t+2 with 1 year lag and matching on missing. It generates the results for row 2 of panel 1 in Table 2.
set.seed(12345)
m2_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 2, forbid.treatment.reversal = FALSE)
mset.m2_rhrv <- m2_rhrv$att
PE.m2_rhrv  <- PanelEstimate(sets = m2_rhrv, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+3 with 1 year lag and matching on missing. It generates the results for row 3 of panel 1 in Table 2.
set.seed(12345)
m3_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 3, forbid.treatment.reversal = FALSE)

mset.m3_rhrv  <- m3_rhrv$att
PE.m3_rhrv   <- PanelEstimate(sets = m3_rhrv, data = rhrv_cont)

summary(PE.m1_rhrv)
summary(PE.m2_rhrv)
summary(PE.m3_rhrv)

##PANEL 2

#This model estimates the effect of sustained targeting on insurgent conscription at time t+1 with 1 year lag and no matching on missing. It generates the results for row 1 of panel 1 in Table 2.
set.seed(12345)
m4_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = FALSE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 1, forbid.treatment.reversal = FALSE)
mset.m4_rhrv <- m4_rhrv$att
PE.m4_rhrv  <- PanelEstimate(sets = m4_rhrv, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+2 with 1 year lag and no matching on missing. It generates the results for row 2 of panel 1 in Table 2.
set.seed(12345)
m5_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = FALSE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +   I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 2, forbid.treatment.reversal = FALSE)
mset.m5_rhrv <- m5_rhrv$att
PE.m5_rhrv  <- PanelEstimate(sets = m5_rhrv, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+3 with 1 year lag and no matching on missing. It generates the results for row 3 of panel 1 in Table 2.
set.seed(12345)
m6_rhrv <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = FALSE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 3, forbid.treatment.reversal = FALSE)

mset.m6_rhrv  <- m6_rhrv$att
PE.m6_rhrv   <- PanelEstimate(sets = m6_rhrv, data = rhrv_cont)

summary(PE.m4_rhrv)
summary(PE.m5_rhrv)
summary(PE.m6_rhrv)

##PANEL 3
set.seed(12345)
m7_rhrv <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                      +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec, 1:2)),
                      size.match = 10, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 1, forbid.treatment.reversal = FALSE)

mset.m7_rhrv <- m7_rhrv$att
PE.m7_rhrv <- PanelEstimate(sets = m7_rhrv, data = rhrv_cont)

set.seed(12345)
m8_rhrv <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                      +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec, 1:2)),
                      size.match = 10, qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 2, forbid.treatment.reversal = FALSE)

mset.m8_rhrv <- m8_rhrv$att
PE.m8_rhrv <- PanelEstimate(sets = m8_rhrv, data = rhrv_cont)

set.seed(12345)
m9_rhrv <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                      +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec, 1:2)),
                      qoi = "att"
                      , outcome.var = "forcedrec",
                      lead = 3, forbid.treatment.reversal = FALSE)

mset.m9_rhrv <- m9_rhrv$att
PE.m9_rhrv <- PanelEstimate(sets = m9_rhrv, data = rhrv_cont)

summary(PE.m7_rhrv)
summary(PE.m8_rhrv)
summary(PE.m9_rhrv)

### Table C3 ----

summary(lm_rhrv <- lm(forcedrec ~ trt, data=rhrv_cont))
summary(lm_rhrv_sc <- lm(forcedrec  ~trt + dyadsfirstyearconflict + intensitynewnum + irregular2 +ideolcomm, data=rhrv_cont))
summary(lm_rhrv_ac <- lm(forcedrec ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + nat_resource_funding, data=rhrv_cont))


###regression clustered on dyad
rhrv_d <- coeftest(lm_rhrv, vcov = vcovCL, cluster = ~dyadid)
rhrv_d_sc <- coeftest(lm_rhrv_sc, vcov = vcovCL, cluster = ~dyadid)
rhrv_d_ac <- coeftest(lm_rhrv_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
rhrv_c <- coeftest(lm_rhrv, vcov = vcovCL, cluster = ~ccode)
rhrv_c_sc <- coeftest(lm_rhrv_sc, vcov = vcovCL, cluster = ~ccode)
rhrv_c_ac <- coeftest(lm_rhrv_ac, vcov = vcovCL, cluster = ~ccode)

stargazer(rhrv_d, rhrv_d_sc, rhrv_d_ac, rhrv_c, rhrv_c_sc, rhrv_c_ac, keep.stat=c("n"), 
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))

##Non IC

### Table C4 ----
set.seed(12345)
m1_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 1, forbid.treatment.reversal = FALSE)
mset.m1_rhrv_nic <- m1_rhrv_nic$att
PE.m1_rhrv_nic  <- PanelEstimate(sets = m1_rhrv_nic, data = rhrv_cont)



#This model estimates the effect of sustained targeting on insurgent conscription at time t+2 with 1 year lag and matching on missing. It generates the results for row 2 of panel 1 in Table 2.
set.seed(12345)
m2_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 2, forbid.treatment.reversal = FALSE)
mset.m2_rhrv_nic <- m2_rhrv_nic$att
PE.m2_rhrv_nic  <- PanelEstimate(sets = m2_rhrv_nic, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+3 with 1 year lag and matching on missing. It generates the results for row 3 of panel 1 in Table 2.
set.seed(12345)
m3_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 3, forbid.treatment.reversal = FALSE)

mset.m3_rhrv_nic  <- m3_rhrv_nic$att
PE.m3_rhrv_nic   <- PanelEstimate(sets = m3_rhrv_nic, data = rhrv_cont)

summary(PE.m1_rhrv_nic)
summary(PE.m2_rhrv_nic)
summary(PE.m3_rhrv_nic)

##PANEL 2

#This model estimates the effect of sustained targeting on insurgent conscription at time t+1 with 1 year lag and no matching on missing. It generates the results for row 1 of panel 1 in Table 2.
set.seed(12345)
m4_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = FALSE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 1, forbid.treatment.reversal = FALSE)
mset.m4_rhrv_nic <- m4_rhrv_nic$att
PE.m4_rhrv_nic  <- PanelEstimate(sets = m4_rhrv_nic, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+2 with 1 year lag and no matching on missing. It generates the results for row 2 of panel 1 in Table 2.
set.seed(12345)
m5_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = FALSE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +   I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 2, forbid.treatment.reversal = FALSE)
mset.m5_rhrv_nic <- m5_rhrv_nic$att
PE.m5_rhrv_nic  <- PanelEstimate(sets = m5_rhrv_nic, data = rhrv_cont)

#This model estimates the effect of sustained targeting on insurgent conscription at time t+3 with 1 year lag and no matching on missing. It generates the results for row 3 of panel 1 in Table 2.
set.seed(12345)
m6_rhrv_nic <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = FALSE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                          +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) +  I(lag(forcedrec_non_ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 3, forbid.treatment.reversal = FALSE)

mset.m6_rhrv_nic  <- m6_rhrv_nic$att
PE.m6_rhrv_nic   <- PanelEstimate(sets = m6_rhrv_nic, data = rhrv_cont)

summary(PE.m4_rhrv_nic)
summary(PE.m5_rhrv_nic)
summary(PE.m6_rhrv_nic)

##PANEL 3
set.seed(12345)
m7_rhrv_nic <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                          +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec_non_ic, 1:2)),
                          size.match = 10, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 1, forbid.treatment.reversal = FALSE)

mset.m7_rhrv_nic <- m7_rhrv_nic$att
PE.m7_rhrv_nic <- PanelEstimate(sets = m7_rhrv_nic, data = rhrv_cont)

set.seed(12345)
m8_rhrv_nic <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                          +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec_non_ic, 1:2)),
                          size.match = 10, qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 2, forbid.treatment.reversal = FALSE)

mset.m8_rhrv_nic <- m8_rhrv_nic$att
PE.m8_rhrv_nic <- PanelEstimate(sets = m8_rhrv_nic, data = rhrv_cont)

set.seed(12345)
m9_rhrv_nic <- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "CBPS.weight",
                          data = rhrv_cont, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + I(lag(fcap, 1:2)) 
                          +  I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(ideolcomm, 1:2)) + I(lag(intensity, 1:2)) + I(lag(nat_resource_funding, 1:2)) +   I(lag(forcedrec_non_ic, 1:2)),
                          qoi = "att"
                          , outcome.var = "forcedrec_non_ic",
                          lead = 3, forbid.treatment.reversal = FALSE)

mset.m9_rhrv_nic <- m9_rhrv_nic$att
PE.m9_rhrv_nic <- PanelEstimate(sets = m9_rhrv_nic, data = rhrv_cont)

summary(PE.m7_rhrv_nic)
summary(PE.m8_rhrv_nic)
summary(PE.m9_rhrv_nic)

### Table C5 ----

summary(lm_nonic <- lm(forcedrec_non_ic ~ trt, data=rhrv_cont))
summary(lm_nonic_sc <- lm(forcedrec_non_ic  ~trt + dyadsfirstyearconflict + intensitynewnum + irregular2 +ideolcomm, data=rhrv_cont))
summary(lm_nonic_ac <- lm(forcedrec_non_ic ~ trt + excluded + exclpopn + dyadsfirstyearconflict + fcap + rsup + irregular2 +ideolcomm +intensitynewnum + nat_resource_funding, data=rhrv_cont))


###regression clustered on dyad
nonic_d <- coeftest(lm_nonic, vcov = vcovCL, cluster = ~dyadid)
nonic_d_sc <- coeftest(lm_nonic_sc, vcov = vcovCL, cluster = ~dyadid)
nonic_d_ac <- coeftest(lm_nonic_ac, vcov = vcovCL, cluster = ~dyadid)


###regression clustered on country

#ic 
nonic_c <- coeftest(lm_nonic, vcov = vcovCL, cluster = ~ccode)
nonic_c_sc <- coeftest(lm_nonic_sc, vcov = vcovCL, cluster = ~ccode)
nonic_c_ac <- coeftest(lm_nonic_ac, vcov = vcovCL, cluster = ~ccode)

stargazer(nonic_d, nonic_d_sc, nonic_d_ac, nonic_c, nonic_c_sc, nonic_c_ac, keep.stat=c("n"), 
          covariate.labels=c("Sustained state targeting", "Excluded", "Excluded population",
                             "Dyads first year", "Fighting capacity", "Rebel support", 
                             "Intensity", "Natural resource funding", 
                             "Irregular", "Communist ideology", "Constant"))
## APPENDIX D: Other forms of RHRV ----


rhrv_cont <- rhrv_cont %>%
  mutate(killings = ifelse(rkillings_s >= 1 | coalesce(rkillings_a, 0) >= 1, 1, 0))

rhrv_cont <- rhrv_cont %>%
  mutate(sexual = ifelse(rsexual_s >= 1 | coalesce(rsexual_a, 0) >= 1, 1, 0))


rhrv_cont <- rhrv_cont %>%
  mutate(displace = ifelse(rdisplace_s >= 1 | coalesce(rdisplace_a, 0) >= 1, 1, 0))

set.seed(12345)
m_k <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "CBPS.weight",
                      data = rhrv_cont, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                      +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(killings, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "killings",
                      lead = 1:3, forbid.treatment.reversal = FALSE)
mset.m_k <- m_k$att
PE.m_k  <- PanelEstimate(sets = m_k, data = rhrv_cont)
summary(PE.m_k)

set.seed(12345)
m_s <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "trt", refinement.method = "CBPS.weight",
                  data = rhrv_cont, match.missing = TRUE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(sexual, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "sexual",
                  lead = 1:3, forbid.treatment.reversal = FALSE)
mset.m_s <- m_s$att
PE.m_s  <- PanelEstimate(sets = m_s, data = rhrv_cont)
summary(PE.m_s)

set.seed(12345)
m_d <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                  treatment = "trt", refinement.method = "CBPS.weight",
                  data = rhrv_cont, match.missing = TRUE,
                  covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) 
                  +  I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(ideolcomm, 1:1)) + I(lag(intensity, 1:1)) + I(lag(nat_resource_funding, 1:1)) + I(lag(displace, 1:1)),
                  size.match = 15, qoi = "att"
                  , outcome.var = "displace",
                  lead = 1:3, forbid.treatment.reversal = FALSE)
mset.m_d <- m_d$att
PE.m_d  <- PanelEstimate(sets = m_d, data = rhrv_cont)
summary(PE.m_d)

### REPLICATION FOR FIGURE D1 ----

png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR/rhrvplot.png")
par(mfrow=c(2,2))
plot(PE.m_k, ylab = "Estimated Effect on Killings",
     xlab = "Time", main = NA, ylim=c(-0.9,0.6))
plot(PE.m_s, ylab = "Estimated Effect on SV",
     xlab = "Time", main = NA, ylim=c(-0.9,0.6))
plot(PE.m_d, ylab = "Estimated Effect on Displacement",
     xlab = "Time", main = NA, ylim=c(-0.9,0.6))
dev.off()

## APPENDIX E: Covariate balance ----

#checking to see which method provides the best balance for main analysis

covpm1_none <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                          treatment = "trt", refinement.method = "none",
                          data = IC, match.missing = TRUE,
                          covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                          size.match = 15, qoi = "att"
                          , outcome.var = "ic",
                          lead = 1:3, forbid.treatment.reversal = FALSE)

pm1_mah <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "mahalanobis",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 1, forbid.treatment.reversal = FALSE)

PE.pm1_mah <- PanelEstimate(sets = pm1_mah, data = IC)

balance_scatter(matched_set_list = list(pm1_mah$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "nat_resource_funding", "ic"))


pm1_psm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "ps.match",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1))+ I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 1, forbid.treatment.reversal = FALSE)

PE.pm1_psm <- PanelEstimate(sets = pm1_psm, data = IC)

balance_scatter(matched_set_list = list(pm1_psm$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "ic"))


pm1_cbps <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                       treatment = "trt", refinement.method = "CBPS.weight",
                       data = IC, match.missing = TRUE,
                       covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1))+ I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                       size.match = 15, qoi = "att"
                       , outcome.var = "ic",
                       lead = 1:3, forbid.treatment.reversal = FALSE)

PE.pm1_cbps <- PanelEstimate(sets = pm1_cbps, data = IC)

png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR")
balance_scatter(matched_set_list = list(pm1_cbps$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ideolcomm", "ic"))
dev.off()

pm1_cbpsm <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                        treatment = "trt", refinement.method = "CBPS.match",
                        data = IC, match.missing = TRUE,
                        covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                        size.match = 15, qoi = "att"
                        , outcome.var = "ic",
                        lead = 1, forbid.treatment.reversal = FALSE)

PE.pm1_cbpsm <- PanelEstimate(sets = pm1_cbpsm, data = IC)
balance_scatter(matched_set_list = list(pm1_cbpsm$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "icpret"))

balance_scatter(matched_set_list = list(pm1_cbpsm$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "icpret"))

#creating covariate balance table for main analysis: IC

#1 year pre-treatment with no refinement 

set.seed(12345)
m1_ref<- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "none",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1))+  I(lag(ic, 1:1)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 1, forbid.treatment.reversal = FALSE)
set.seed(12345)


#covariate balance for 1 year pre-treatment, before and after refinement
get_covariate_balance(m1_noref$att, IC, covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ic"), plot = FALSE)
get_covariate_balance(m1$att, IC, covariates =  c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum","nat_resource_funding", "ic"), plot = FALSE)


#2 year pre-treatment with no refinement 

set.seed(12345)
m7_noref<- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                      treatment = "trt", refinement.method = "none",
                      data = IC, match.missing = TRUE,
                      covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + + I(lag(fcap, 1:2)) +  I(lag(ideolcomm, 1:2)) + I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(nat_resource_funding, 1:1))+ I(lag(ic, 1:2)),
                      size.match = 15, qoi = "att"
                      , outcome.var = "ic",
                      lead = 1, forbid.treatment.reversal = FALSE)

#covariate balance for 2 year pre-treatment, before and after refinement
get_covariate_balance(m7_noref$att, IC, covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ic"), plot = FALSE)
get_covariate_balance(m7$att, IC, covariates =  c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ic"), plot = FALSE)



#creating covariate balance table for models including rebel goals

#1 year pre-treatment with no refinement 
set.seed(12345)
PM.sec1_noref<- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                           treatment = "trt", refinement.method = "none",
                           data = IC, match.missing = TRUE,
                           covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(ic, 1:1)) + I(lag(nat_resource_funding, 1:1))+ I(lag(sec, 1:1)),
                           size.match = 15, qoi = "att"
                           , outcome.var = "ic",
                           lead = 1, forbid.treatment.reversal = FALSE)


#covariate balance for 1 year pre-treatment, before and after refinement
get_covariate_balance(PM.sec1_noref$att, IC, covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "ic", "nat_resource_funding", "sec"), plot = FALSE)
get_covariate_balance(PM.sec1$att, IC, covariates =  c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "ic", "nat_resource_funding", "sec"), plot = FALSE)


#2 year pre-treatment with no refinement 

set.seed(12345)
PM.sec2_1_noref<- PanelMatch(lag = 2, time.id = "year", unit.id = "dyadid",
                             treatment = "trt", refinement.method = "none",
                             data = IC, match.missing = TRUE,
                             covs.formula = ~  I(lag(excluded, 1:2)) + I(lag(exclpopn, 1:2)) + I(lag(dyadsfirstyearconflict, 1:2)) + + I(lag(fcap, 1:2)) +  I(lag(ideolcomm, 1:2)) + I(lag(rsup, 1:2)) + I(lag(irregular2, 1:2)) + I(lag(intensitynewnum, 1:2)) + I(lag(ic, 1:2)) + I(lag(nat_resource_funding, 1:1))+  I(lag(sec, 1:2)),
                             size.match = 15, qoi = "att"
                             , outcome.var = "ic",
                             lead = 1, forbid.treatment.reversal = FALSE)

#covariate balance for 2 year pre-treatment, before and after refinement
get_covariate_balance(PM.sec2_1_noref$att, IC, covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "ic", "nat_resource_funding", "sec"), plot = FALSE)
get_covariate_balance(PM.sec2_1$att, IC, covariates =  c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "ic", "nat_resource_funding", "sec"), plot = FALSE)


#creating covariate balance table for models with alternate treatment threshold

# 1 year pre-treatment with no refinement 
set.seed(12345)
am1_noref<- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid",
                       treatment = "alttrt", refinement.method = "none",
                       data = IC, match.missing = TRUE,
                       covs.formula = ~  I(lag(excluded, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + + I(lag(fcap, 1:1)) +  I(lag(ideolcomm, 1:1)) + I(lag(rsup, 1:1)) + I(lag(irregular2, 1:1)) + I(lag(intensitynewnum, 1:1)) + I(lag(nat_resource_funding, 1:1))+ I(lag(ic, 1:1)),
                       size.match = 15, qoi = "att"
                       , outcome.var = "ic",
                       lead = 1, forbid.treatment.reversal = FALSE)


#covariate balance for 1 year pre-treatment, before and after refinement
get_covariate_balance(am1_noref$att, IC, covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ic"), plot = FALSE)
get_covariate_balance(am1$att, IC, covariates =  c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "ideolcomm", "rsup", "irregular2", "intensitynewnum", "nat_resource_funding", "ic"), plot = FALSE)



#scatter plots
# Main models
png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR/cbpsbalancemain.png")
balance_scatter(matched_set_list = list(m1$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "nat_resource_funding", "ic"))

dev.off()



#Including rebel goals

png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR/cbpsbalancesec.png")
balance_scatter(matched_set_list = list(PM.sec1$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "ic", "nat_resource_funding", "sec"))
dev.off()


#alternate treatment

png("/Users/emilymyers/Library/CloudStorage/Box-Box/Research Projects/Insurgent Conscription Project/JCR Submission/Data Analysis/Replication data_JCR/cbpsbalancealttrt.png")
balance_scatter(matched_set_list = list(am1$att),
                data = IC,
                covariates = c("excluded", "exclpopn", "dyadsfirstyearconflict", "fcap", "rsup", "irregular2", "intensitynewnum", "ideolcomm", "nat_resource_funding", "ic"))
dev.off()







